What if we wanted to look at population in 20 years given an initial condition
Two options
explicit solution to differential equation is known; e.g. you can integrate both sides of the equation! Not always possible but lets look at a case where it is possible
must be solved by iteration; this is what we do when we can’t integrate both sides
## function (T, P0, r, K)
## {
## P = P0 * exp(r * T)
## if (P > K) {
## P = K
## }
## return(P)
## }
# gives population after any time given an initial population
# 20 rabbits, growth rate of 0.01 how many in 30 years
exppop(T=30, P0=20, r=0.01, K=1000)## [1] 26.99718
# if we want to see how population evolves over time - generate a time series by running our model for each point in time
initialrabbits = 20
years = seq(from=1, to=100, by=2)
Ptime = years %>% map_dbl(~exppop( P0=initialrabbits, r=0.01, K=1000, T=.x))
# keep track of what times we ran
Ptime = data.frame(P=Ptime, years=years)
ggplot(Ptime, aes(years,P))+geom_point()+labs(x="years",y="Rabbit Population")Population models can be discrete (rather than continuous)
So we could implement them as difference equations and iterate
source("../R/discrete_logistic_pop.R")
# notice how a for loop is used to iterate
# how many rabbits after 50 years given a growth of 0.1
# starting with 1 rabbit - but a carrying capcity of 500
discrete_logistic_pop## function (P0, r, K, T = 10)
## {
## pop = P0
## for (i in 1:T) {
## pop = pop + r * pop
## pop = ifelse(pop > K, K, pop)
## }
## return(pop)
## }
## [1] 11.4674
## [1] 12.18249
## [1] 12.18249
## [1] 11.4674
# why are they different
# look at trajectories
growth_result = data.frame(time=seq(from=1,to=100))
growth_result$Panalytic = growth_result$time %>% map_dbl(~exppop( P0=1,r=0.05, K=200,T=.x ))
growth_result$Pdiscrete = growth_result$time %>% map_dbl(~discrete_logistic_pop( P0=1,r=0.05, K=200,T=.x ))
tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()Using a solver….when you can’t do the integration by hand :)
For example, if you added a carrying capacity as threshold where it stops growing
In this case we integrate by iteration - approximates analytic integration
Implement the differential equation as a function that
*inputs time, the variable(s) and a parameter list
(it needs time even though you don’t use it)
My convention: name derivative functions starting with d to remind myself that they are computing a derivative
Only works for Ordinary Differential Equations - single independent variable (in our case time)
Partial differential equations - more than 1 independent variable (e.g x and y if changing in space)
R has a solver called ODE for solving ordinary differential equations from package desolve
## function (time, P, r)
## {
## dexpop = r * P
## return(list(dexpop))
## }
library(deSolve)
initialrabbits = 20
years = seq(from=1, to=100, by=2)
# run the solver
Ptime = ode(y=initialrabbits, times=years, func=dexppop,parms=c(0.01))
head(Ptime)## time 1
## [1,] 1 20.00000
## [2,] 3 20.40404
## [3,] 5 20.81623
## [4,] 7 21.23675
## [5,] 9 21.66576
## [6,] 11 22.10344
colnames(Ptime)=c("year","P")
# notice that there are additional pieces of information year, including the method used for integration
attributes(Ptime)## $dim
## [1] 50 2
##
## $dimnames
## $dimnames[[1]]
## NULL
##
## $dimnames[[2]]
## [1] "year" "P"
##
##
## $istate
## [1] 2 52 105 NA 6 6 0 36 21 NA NA NA NA 0 1 1 NA NA NA
## [20] NA NA
##
## $rstate
## [1] 2.00000 2.00000 99.98839 0.00000 1.00000
##
## $lengthvar
## [1] 1
##
## $class
## [1] "deSolve" "matrix"
##
## $type
## [1] "lsoda"
# this also means you need to extract just the data frame for plotting
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")# this also works (of course function can be by order)
Ptime=ode(initialrabbits, years, dexppop,0.03)
colnames(Ptime)=c("year","P")
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")You can play a bit with changing your function to something that you can’t integrate “by hand”
BUT we might want more parameters
to work with ODE, parameters must all be input as a single list; similar to how we return multiple outputs from a function
see example below..lets add a carrying capacity
## function (time, P, parms)
## {
## dexpop = parms$r * P
## dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
## return(list(dexpop))
## }
# create parameter list
initalrabbits=2
newparms = list(r=0.03, carry_capacity=300)
#apply solver
results=ode(initialrabbits, years, dexppop_play,newparms)
head(results)## time 1
## [1,] 1 20.00000
## [2,] 3 21.23675
## [3,] 5 22.54993
## [4,] 7 23.94435
## [5,] 9 25.42499
## [6,] 11 26.99719
# add more meaningful names
colnames(results)=c("year","P")
# plot
ggplot(as.data.frame(results), aes(year,P))+geom_point()+labs(y="Population", "years")# try again with different parameters
alternativeparms = list(r=0.04, carry_capacity=500)
results2=ode(initialrabbits, years, dexppop_play,alternativeparms)## DLSODA- At current T (=R1), MXSTEP (=I1) steps
## taken on this call before reaching TOUT
## In above message, I1 = 5000
##
## In above message, R1 = 81.5108
##
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
## time 1
## [1,] 1 20.00000
## [2,] 3 21.66575
## [3,] 5 23.47022
## [4,] 7 25.42499
## [5,] 9 27.54257
## [6,] 11 29.83650
colnames(results2)=c("year","P_parm2")
# plot
ggplot(as.data.frame(results2), aes(year,P_parm2))+geom_point()+labs(y="Population", "years")# compare by combining into a single data frame
both = inner_join(as.data.frame(results), as.data.frame(results2))## Joining, by = "year"
both_p = both %>% gather(key="model", value="P", -year)
ggplot(both_p, aes(year,P, col=model))+geom_point()+labs(y="Population", "years")Remember we have 3 ways now to calculate population
analytical solution - based on integration (exppop.R) BEST
using an ode solver for numerical approximation (exppop_play.R)
numerical integration using in discrete steps (discrete_logistic_pop.R)
Finally look at continuous derivative using ODE solve Needs initial conditions differential equation *parameters
## function (time, P, parms)
## {
## dexpop = parms$r * P
## dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
## return(list(dexpop))
## }
# set up using the same parameters
pcompare = list(r=r,carry_capacity=K)
# now run our ODE solver
result = ode(P0, growth_result$time, dexppop_play, pcompare)
head(result)## time 1
## [1,] 1 1.000000
## [2,] 2 1.051273
## [3,] 3 1.105172
## [4,] 4 1.161835
## [5,] 5 1.221404
## [6,] 6 1.284027
# we already have time - so just extract population
growth_result$Pdifferential = result[,2]
# compare all 3 approaches
tmp = growth_result %>% pivot_longer(cols=-time,names_to="Ptype",values_to="P")
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()All differential and difference equations are approximations The analytical solution is exact
Notice that differential equations is a bit more accurate!
Diffusion can be implemented as a partial differential equation
Complicated to solve - but there are tool in R for specific types of partial differential equations
More info on differential equations in R
Diffusionn would require partial derivatives - time and space! it gets much more tricky …beyond this course
But we can appoximate diffusion with a difference equation - and iterative to get an estimate of how diffusion works
Example of Diffusion - difference equation implementation to see what some issues can be